Computing stationary free-surface shapes in microfluidics 
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A finite-element algoritlim for computing free-surface flows driven by arbitrary body forces is 
presented. The algorithm is primarily designed for the microfluidic parameter range where (i) the 
Reynolds number is small and (ii) force-driven pressure and flow fields compete with the surface 
tension for the shape of a stationary free surface. The free surface shape is represented by the 
boundaries of finite elements that move according to the stress applied by the adjacent fiuid. Ad- 
ditionally, the surface tends to minimize its free energy and by that adapts its curvature to balance 
the normal stress at the surface. The numerical approach consists of the iteration of two alternating 
steps: The solution of a fluidic problem in a prescribed domain with slip boundary conditions at the 
free surface and a consecutive update of the domain driven by the previously determined pressure 
and velocity fields. For a Stokes problem the first step is linear, whereas the second step involves the 
nonlinear free-surface boundary condition. This algorithm is justified both by physical and mathe- 
matical arguments. It is tested in two dimensions for two cases that can be solved analytically. The 
magnitude of the errors is discussed in dependence on the approximation order of the finite elements 
and on a step-width parameter of the algorithm. Moreover, the algorithm is shown to be robust 
in the sense that convergence is reached also from initial forms that strongly deviate from the final 
shape. The presented algorithm does not require a remeshing of the used grid at the boundary. 
This advantage is achieved by a built-in mechanism that causes a smooth change from the behavior 
of a free surface to that of a rubber band if the boundary mesh becomes irregular. As a side effect, 
the element sides building up the free surface in two dimensions all approach equal lengths. The 
presented variational derivation of the boundary condition corroborates the numerical finding that 
a second-order approximation of the velocity also necessitates a second-order approximation for the 
free surface discretization. 

PACS numbers: 47.55.Ca, 47.61. Jd, 47.11. Fg, 47.10.A-, 47.15.G- 



I. INTRODUCTION 

In the past decade the development of so-called 
"labs-on-a-chip" has led to an increased interest in 
microfluidics,^"^ i.e. in the field of hydrodynamics with 
characteristic length scales of less than a millimeter. 
These flows are characterized by small Reynolds numbers 
and consequently governed by the Stokes equations. In 
the case of prescribed fluid domains with no-slip bound- 
ary conditions standard numerical methods exist for com- 
puting their solutions. ^'^ 

Recently, various experimental techniques®'^ have been 
developed to induce and control flows in fluids which sit 
on a substrate without being conflned by lateral and cov- 
ering walls. ^"'^^ In the experiments the fluid is kept to- 
gether by its surface tensions both at the substrate and 
at the fluid-air interface. The stationary form that is as- 
sumed by the fluid-air interface is not given a priori. It 
results from an interplay of the internal streaming pat- 
tern, the internal pressure distribution and the surface 
tension. On the other hand, the form of the interface 
acts back on the flow. This mutual interaction of form 
and flow renders free boundary value problems fascinat- 
ing but difficult. The relative importance of viscous fiow 
and pressure, each compared to the influence of the sur- 
face tension, can be quantifled by two dimensionless num- 
bers, the capillary number and a generalized Bond num- 
ber^ respectively. 

In the present work we consider a small water droplet 



(around 50 nl or less). The droplet sits on a flat substrate 
and is mechanically agitated by a body force. Inside the 
droplet this body force then causes stationary pressure 
and flow fields which can lead to a significant deforma- 
tion for the free surface. Sufficiently strong body forces 
may lead to the motion of the entire droplet, but this sit- 
uation will not be considered here. The values of all, the 
Reynolds number, the capillary number, and the Bond 
number are assumed to range from zero up to unity. This 
corresponds to experimentally relevant situations. ^"'^^'^^ 

Several numerical approaches for determining free sur- 
face shapes have been proposed in the past. The suitabil- 
ity of the approaches depends on the size of the system, 
typical velocities, and the material properties, as well as 
on the resulting deformation of the fluid domain. They 
can roughly be classifled into two groups: Either a flxed 
grid and a function describing the position of the free 
surface is used, or the computational mesh is moved to- 
gether with the fluid domain, yielding a sharp surface 
representation by elements' boundaries. 

An established method of the first kind is the con- 
tinuum method proposed by Brackbill et. al}^ They cir- 
cumvented the discretization of the normal-stress bound- 
ary condition by introducing a body force density that 
is concentrated near the free surface. This force den- 
sity accounts for the effect of surface tension. We have 
tested this method, which is implemented in the commer- 
cially available fluid-dynamics program FLUENT using 
a volume-of-fluid discretization. For a macroscopic sys- 
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tern this method worked fine. The method, however, fails 
if the system is scaled down to the microfluidic param- 
eter regime. In a simple test example we found that 
approximation errors of the free-surface boundary con- 
dition contributed to the force balance in the Navier- 
Stokes equations and were amplified in an uncontrolled 
manner. This typically gave rise to a spurious veloc- 
ity field. It even occurred when we started the iteration 
with the known solution. Problems with this method 
have also been reported by Renardy & Renardy^^ and 
by Popinet & Zaleski.^^ Lafaurie et. al}"^ find the spuri- 
ous velocities to be of the order surface-tension/viscosity 
which is the dominant velocity scale for microfluidic sys- 
tems. Thus, the existing continuum method appears to 
be inappropriate for the microfluidic parameter regime. 

Another approach of the flrst kind has recently been 
proposed by Smolianski.^* He uses finite elements and 
a level-set description for the free surface and calculates 
curvatures by derivatives of the distance-function. He 
still encounters spurious velocity fields proportional to 
the ratio surface- tension/viscosity. 

Methods of the second kind, representing the free sur- 
face by a sharp interface arc expected to work better in 
the microfiuidic parameter regime. Algorithms in this 
class are often referred to as "moving mesh" or "ALE" 
methods and generally require more involved techniques, 
keeping the computational mesh feasible and not too dis- 
torted. 

A technique of the second kind that has successfully 
been employed for tension-dominated free-surface prob- 
lems is the boundary-element method. The dimen- 
sionality of the equations is reduced to the dimension- 
ality of the surface which provides the basis for an effi- 
cient implementation. Unfortunately, this reduction can 
only be performed for Stokes equations with conserva- 
tive body forces, which can be absorbed into the pres- 
sure term. In the present investigation we allow for non- 
conservative body forces which are of particular experi- 
mental relevance. ^'^'^ 

Pioneering works for the finite-element implementation 
of the full free-surface problem were published by Scriven 
and coworkers. ^^'^^ They used spines to parameterize the 
movements of the compiitational mesh in coating flow 
and implemented Newton's method for a Galerkin ap- 
proximation scheme. This work was later continued un- 
der the designation "total linearization method" by Cu- 
velier and coworkers. ^'^'^^ Their description requires a 
height function for the free-surface position, which makes 
it necessary to use well-adapted coordinate systems like 
polar cylindrical or spherical ones. It must be known in 
advance if a free surface will overhang. 

In the present paper we extend the works of Scriven 
and Cuvelier to arbitrary surface geometries. In our 
description, the parameterization of the free surface is 
given directly by the finite-elements' boundary parame- 
terization. Thus, neither spines nor a height function are 
needed. To properly account for intrinsic curvatures of 
the free surface, all equations are formulated in a fully 



covariant form that allows for all differential-geometric 
properties of the surface. An excellent reference for this 
formulation are the works of Aris^^ and Scriven^^ where 
the fluidic flow inside a curved free surface is described. 

Recently, algorithms have been published that describe 
fully time-dependent free-surface flows, even in three 
dimensions. ^^"^^ In these works the free surface is moved 
mainly due to the kinematic boundary condition, i.e. it 
is advected passively. Concerning convergence, there has 
been a controversy if the kinematic or rather the nor- 
mal stress boundary condition should be used to move 
the free surface. This issue was resolved by Silliman and 
Scriven who state that for capillary numbers below unity, 
the normal stress iteration converges well while a kine- 
matic iteration eventually fails. ^-'^ In addition, when the 
kinematic boundary condition is used for updating the 
free surface, the balance of normal stress that carries the 
effects of surface tension is not strictly imposed. It is 
used when implementing the weak form of the Navier- 
Stokes equations: In this context an integration by parts 
yields an integral of the normal stress over the free sur- 
face, which is then replaced by the corresponding surface- 
integral of the tension forces. Similar techniques are com- 
monly used for problems with outflow boundary condi- 
tions or for Poisson's equation with Neumann boundary 
conditions. The correctness of the technique has been 
justified for the outfiow problem by Renardy. However, 
it is not evident if it also works in the case where the 
surface-tension terms dominate the whole problem. The 
question remains open, in which sense the boundary con- 
dition is satisfied. Therefore, we found it necessary in 
our examples to visualize the terms involved in the free- 
surface boundary condition, thus proving that they are 
correctly balanced. 

An important result of a variational description of the 
tension terms is an improvement of the Newton algorithm 
controlling possible mesh distortions at the free surface. 
Many algorithms implementing the weak form of the cap- 
illary boundary condition encounter intrinsic instabilities 
of the boundary mesh when significant changes of the free 
surface take place. For the program surface evolver^^ this 
manifests itself in shrinking and growing surface facets. 
Similar effects have been observed by Brinkmann'^^ and 
Bansch.^^ Our formulation of the capillary free surface 
is such that the free surface smoothly changes to the be- 
havior of a rubber band when the boundary mesh be- 
comes distorted. This leads to an automatic regulariza- 
tion of the mesh without the need of explicit remeshing 
or smoothing. 

In Section II the mathematical formulation of the prob- 
lem is presented in terms of differential equations, to- 
gether with the boundary conditions and the relevant 
parameter regime. In Sec. Ill we then re-establish the 
bulk equations and their boundary conditions by vari- 
ational techniques. For the free-surface we introduce a 
differential-geometric notation that allows us to write the 
boundary condition in a weak form. Up to this point a 
continuous description is used. Section IV introduces the 
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discretization of the problem by the computational mesh. 
The formulation of the tension forces as the concurrent 
minimization of the free-surface area of single finite ele- 
ments is a necessary requirement for the mentioned auto- 
matic regularization mechanism. Sec. V provides a short 
summary of the whole algorithm. In Sec. VI we present 
examples that show the accuracy of the algorithm and 
two further examples for different values of the capillary 
and Bond numbers. Mathematical and algorithmic de- 
tails are deferred into the Appendices A-D. 



II. STATEMENT OF THE PROBLEM 

Throughout the paper, we shall write all equations in 
tensor notation for arbitrary curvilinear coordinate sys- 
tems. This will considerably simplify the differential geo- 
metric notation in the following sections. For the formu- 
lation of the full Navier-Stokes equations in curvilinear 
coordinates we refer to Aris.^^ A repeated index that oc- 
curs in CO- and contravariant positions is summed over, 
indices that are preceded by a comma denote covariant 
derivatives, and gij is the metric tensor of the underlying 
coordinate system. 



A. The basic equations 

We study incompressible and stationary flows, that are 
characterized by a small Reynolds number Re = pxv/rj. 
Here, p is the density of the fluid, r/ its viscosity, and 
X and V denote typical magnitudes of length and veloc- 
ity. Under these conditions the pressure field p and the 
velocity field with components satisfy the Stokes equa- 
tions^^ 



<■ = 0, 



(1) 



where aij is the fluidic stress tensor, /* an external body 
force causing non-trivial streaming and pressure patterns 
within a domain V. It can be split into its conserva- 
tive part /* which can be displayed as the gradient of 
a potential and its non-conservative part f^^ with van- 
ishing divergence. The domain V may be bounded by 
rigid walls and by free surfaces, as e.g. a droplet sitting 
on a substrate. Equations (1) and (2) then are subject 
to boundary conditions at different parts of the bound- 
ary dV: First, the flow has to meet the kinematic bound- 
ary condition, requiring that the normal projection of a 
stationary velocity field vanishes at the boundary, i.e., 



ViN' = 0. 



(3) 



At immobile sticky walls we use the no-slip boundary 
condition, according to which the velocity vanishes also 
in the tangential directions of the boundary, implying 



Here, denotes the ith component of the tangential 
vector Tq (a ~ 1,2 for a two-dimensional surface). The 
remaining boundary is a free surface that dynamically 
adjusts its position such that the stress balance holds. 



cr'-' Nj = jkN^ on free surfaces, 



(5) 



with the surface tension 7 and the curvature k. Note 
that we have omitted the gradient of the surface tension 
and thus exclude Marangoni eflFects. This simplifies the 
following calculations but does not present a principal 
restriction of our description. 

Equations (l)-(5) have been simplified by assuming 
that the fiuid in domain V is surrounded by a medium 
of much smaller viscosity, which is the case e.g. for a 
water— air interface at room temperature. Therefore, the 
surrounding's viscous stress contribution does not show 
up in the balance Eq. (5). We further assume that the 
ambient pressure pq is homogeneous. Since the pressure 
is determined by the Stokes equations only up to a con- 
stant, we can split it into a part pi with vanishing average 
and use the ambient pressure po as an offset parameter 
which enters only in the normal stress balance (5), 



p(x) =po+pi (x) with 



/ pi dV = 0. 
Jv 



(6) 



B. The parameter regime 



By transforming both, the bulk equation (2) and the 
free boundary condition (5) into dimensionless form em- 
ploying viscosity-scaling, one observes that the system 
may be characterized by two relevant ratios of forces, 
given by the dimensionless numbers 



Bo 
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and Ca = 



7 " 



(7) 



Here, x,v and fc denote typical magnitudes of length, 
velocity and the conservative part of the force density, 
respectively. Bo is a generalization of the Bond number 
which is usually defined in terms of gravitational forces 
only. The capillary number Ca measures the viscous con- 
tribution to the surface deformation. In a system with 
static boundaries and vanishing Reynolds number we can 
express the velocity scale by the typical magnitude fnc 
of the non-conservative part of the driving force, namely 
V = f^c/rj. This yields an alternative definition of the 
capillary number similar to that of the Bond number, 



Ca = 



fnc^ 

7 



(8) 



ViTl = at the walls. 



(4) 



These two numbers reflect the very different effects of 
the conservative and the non-conservative parts of the 
driving. In this sense, Ca provides also a measure for the 
spatial changes of the velocity field. For small Ca the 
flow is slow and changes smoothly, whereas for large Ca 
it may exhibit drastic gradients. 
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FIG. 1: A sketch of the coordinate system on a two- 
dimensional surface A, embedded into the three-dimensional 
space. The surface coordinates v are mapped from the refer- 
ence domain E (left) onto the surface A (right) via the par- 
rameterization vector t(i/). 



We propose our numerical scheme for the parameter 
regime where both, Ca, and Bo are of order unity or less. 
Thus, pressure gradients and viscous forces can deform 
the free surface significantly. The surface tension is large 
enough, however, in order to keep the whole fluid domain 
together, a pinch-off cannot occur. The viscosity renders 
the velocity field smooth over the whole fluid domain 
and prevents the existence of boundary layers. Because 
we consider only stationary flows in stationary domains 
according to Eqs. (2) and (3), and because we have set the 
substrate's velocity to zero in Eq. (4), we always obtain 
pinned contact-lines. A rolling or slipping droplet would 
raise additional challenges regarding the stress near the 
contact-line that are beyond the scope of this paper. 



III. CONTINUOUS DESCRIPTION OF THE 
PROBLEM 

In order to clarify the numerical treatment of the free- 
surface boundary condition we flrst explore the physical 
origins of the balanced forces. Wc will then express each 
of them by the flrst variation of a functional. For New- 
ton's method it will be necessary to calculate also the 
second variation. 



Physical aspects of the free-surface boundary 
condition: first variations 



The surface tension term "fKNi in the boundary con- 
dition (5) arises from the fact that an extended inter- 
face between two different phases "costs" free energy. ^'^ 
To flnd the optimal configuration the surface is continu- 
ally probing positions in its vicinity in order to minimize 
its free energy. For the case of an applied conservative 
force fi = — the system is static (t;* = 0), and the 
free-surface boundary condition is equivalent to a mini- 
mization of a free energy expression. This calculation is 
performed in Appendix A. 



Due to its thermodynamic origin, the surface tension 
term results from a first variation of a functional. This 
carries over also to the dynamic case (w* ^ 0) in which the 
boundary condition (5) must hold at any instant of time. 
The contribution of the free surface A to the system's free 
energy is given by the integral of the surface tension 7 
over A, where dA denotes the infinitesimal surface area, 



F = 



jdA. 



(9) 



Any smooth surface in a D-dimensional space may be 
parameterized by D—1 surface coordinates v" (a = 
1) which determine the coordinates t*(i^") of 
points in D-dimensional space on the surface. Both, sur- 
face and space coordinates are illustrated in Fig. 1. In 
our numerical studies wc restrict ourselves to D = 2. The 
general framework, however, remains valid also for D = 
3. The surface coordinates j^" are taken from the pa- 
rameter set E C R^~^. With v running through E the 
whole free surface A is covered. 



A = {e^i/i^)\ueE}. 



(10) 
(11) 



Here, e(j) is the ith base vector in space. Throughout 
the paper we will always use Greek symbols for surface 
indices and Latin ones for space indices. The connection 
between surface and space coordinates is conveniently de- 
scribed by the surface-derivatives of the parameterization 
functions (cf. Aris,^^ p. 215), i. e., 



(12) 



Understood as a contravariant £)-dimensional space- 
vector, t*^ represent the components of the ath tangent 
vector Ta to the surface. At the same time, is a 
covariant surface- vector. We can now construct the com- 
ponents of the surface's metric tensor the scalar 
products of these tangent vectors, namely 

aa/3 = 9ij t]^ t^^0 and a = det{aa0). (13) 

The metric tensor Uaff, its determinant a and its in- 
verse a"^ are nonlinear functions of the tangential vector 
components t^^, see in Appendix B. The normal vector 
the normalized cross-product of two tangent vectors, 



/3' 



(14) 



and the curvature k is given as the trace of the tensor ba/3 
of the second fundamental form of the surface, 



K = a"^bai3 with 



(15) 
(16) 



Using the parameterization (10) of the free surface, we 
demonstrate in Appendix B that the change of the free 
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energy contribution F with respect to a variation of the 
surface vectors is given by 

^F[^t] = I^Jt^p 9ij dA. (17) 

By an integration by parts this expression can be cast 
into a form containing the curvature term of the free- 
surface boundary condition (5), i.e., 

6F[6t] = - j -fKNiSf dA (18) 

(sec Appendix C also for the case of varying surface ten- 
sion). In this way, the curvature term n in Eq. (18) 
that contains second spatial derivatives is replaced by 
a product of two terms, each containing a first deriva- 
tive in Eq. (17). Especially for numerical applications it 
is much more favorable to work only with first deriva- 
tives. This trick has been used in the literature in dif- 
ferent contexts. ^"'^^'^""^'^^'^^ Seen from a physical perspec- 
tive, version (17) of the equation is the more natural one. 
Here, one directly deduces that forces pulling along the 
tangential direction attempt to minimize the facet area 
of a mesh's boundary. On the basis of single finite ele- 
ments this perspective will be used below for stabilizing 
the computational mesh. 

The left-hand side of Eq. (5), (TijN-\ is the normal 
fluidic stress at the boundary. We now recapitulate how 
this term can be understood as the result of a variational 
principle. In a stationary system with rigid immobile 
boundaries there is a balance between the power output 
due to viscous dissipation and the power input due to 
external driving. The total power output of the fluid, 
which we will denote with 

P= j VdV= I {a'^Vij- fvi)dV, (19) 
Jv Jv 

vanishes. Additionally, for a prescribed domain V , the 
Stokes equations yield those velocity and pressure fields 
that render the local power extremal (see Finlayson,^^ 
p. 271). Vice versa, from the vanishing first variation 
of P, 

5P[5p\^ f ^SpdV^- I v\6pdV (20) 
Jv cip Jv 

= I {-f6vi + a'^6vij)dV (22) 
Jv 

= -/(/' + '^'i)Svi dV+l a'^NjSvi dA, (23) 
Jv JdV 

the Stokes equations follow by setting the bulk contri- 
butions to zero. The boundary integral in the last row 
of Eq. (23) provides the fluidic stress in the free-surface 
boundary condition. 



At this point we see that the two terms in the stress 
balance Eq. (5) have different physical origins. The sur- 
face tension is of thermodynamic (or rather of "thermo- 
static") nature while the fluidic stress stems from dy- 
namic considerations. The first results from minimizing 
a free energy, while the second stems from minimizing a 
power. Formally, this is expressed by the different varia- 
tions Sv^ and Sf^ in the expressions aijN^Sv'^ in Eq. (23) 
and K"/NiSt^ in Eq. (18). Already for dimensionality rea- 
sons they cannot be equal, neither can the functionals 
F and P be directly combined into one single variational 
principle. 

From an algorithmic point of view one has to make 
a choice here: to approximate the free-surface boundary 
condition using either the Sv^ or the Sf^ as test-functions. 
In our Galerkin implementation of the problem we will 
use the ansatz functions as test functions. Therefore, in 
order to acquire a consistent numerical algorithm we have 
to approximate both, velocity and the geometry param- 
eterization with finite elements of the very same order. 
This is the first central statement of the present work. 

It was stated by Bansch^'' (p. 42, cf. also citations 
49 and 50 therein) that a second-order approximation of 
the surface parameterization yields a "good discrete cur- 
vature" , whereas a first-order one does not. The same can 
be seen below in Fig. 3. We are now able to substantiate 
his numerical observation with the underlying physical 
mechanism. The argument is similar to that for the cele- 
brated Ladyzhenskaya-Babuska-Brezzi requirement that 
velocity gradients have to be approximated by the same 
order as the pressure. From a physical perspective this 
is not astonishing, because both are components of the 
same stress tensor. 



B. Splitting the problem into two numerical 
systems 

For free boundaries a twofold problem must be solved: 

(i) The unknown fluid domain V is to be determined and 

(ii) the Stokes Eqs. (1) and (2) are to be solved within V, 
using the boundary conditions (3)-(5). The latter them- 
selves depend on the shape of V via the normal vector 
at the boundary. Both parts of this problem cannot be 
processed independently. 

In principle, there exist two options to deal with this 
combined problem. A first one is to implement a sin- 
gle numerical system for both, the flow variables p and 

together with the geometry variables t^. We will not 
follow this direction but rather consecutively solve two 
smaller systems, one for the flow variables, depending 
on the current domain V, and a second one for the pa- 
rameterization of the boundary. We have chosen this 
approach because the problem is linear in the flow vari- 
ables and highly nonlinear in the geometry variables t^. 
The nonlinearity is due to the appearance of the inverse 
surface metric a^"^ in Eq. (17). Thus, solving the Stokes 
equations in the fluidic system, how we will call it, will 
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be a standard problem, while the nonlinear search for 
the correct boundary shape will be done in the geometry 
system. Both systems are solved consecutively: 

1. Choose an initial domain V . 

2. Until convergence repeat the following steps: 

(a) Solve the fluidic system within the domain V . 

(b) Solve the geometry system using fixed values 
for the pressure and velocity variables. This 
results in an updated domain V. 

In three-dimensional space the equations (3) (5) pose 
four boundary conditions. They are one too many for 
the linear fluidic system to be fully determined. One 
boundary condition is thus used for updating the param- 
eterization of the free surface. ^"^ The main challenge is the 
proper assignment of specific boundary conditions to the 
two systems in order to make them solvable, uniquely de- 
termined, and robust. It is clear that the no-slip bound- 
ary condition (4) at sticky walls applies only to the fluidic 
system. The free-surface boundary condition yet needs 
further consideration. 

Here, again, a physical argument helps to choose the 
proper boundary condition. It is either the stress by 
the fluid or its velocity that is moving the free surface. 
Accordingly, either the normal stress balance (5) or the 
kinematic boundary condition (3) can be used by the ge- 
ometric system to update the surface (see the discussion 
by Saito & Scriven^^ and our remarks in the introduc- 
tion). We choose our approach according to the following 
principle: The fluidic system should be well deflned as a 
stationary system even if the boundary is fixed and is 
not part of the problem. If then the kinematic boundary 
condition were not imposed on the stationary flow, the 
velocity fleld would pass through the free surface which 
is also stationary. This excludes surface-updates by the 
kinematic boundary condition. 

Until the correct boundary shape has been found, it is, 
in principle, possible that the surrounding flow forces the 
free surface into an arbitrary direction. By its very na- 
ture, however, the tension force stays always normal on 
the free surface. Only normal forces can be compensated 
by a free surface. As a necessary condition the tangential 
projection of the normal stress has to vanish. Whenever 
tangential components emerge during the run of an al- 
gorithm, the result will be a numerical artefact. In the 
proposed scheme with two separated systems it is the flu- 
idic system which must ensure the tangential components 
of the free boundary condition, i.e., 

{vi,j + Vj,i)NH^^^ = foraUa. (24) 

Here, the surface tension 7 has been set constant along 
the surface. For the velocity variables this constitutes 
a perfect slip boundary condition, which is similar to a 
Neumann boundary condition. We thus find the fluidic 
system to be fully determined and physically well defined 



even for fixed boundaries by the conditions (3), (4) and 
(24). 

The geometry system is then responsible for the re- 
maining normal component of the stress balance (5), 

-p + rjivij + Vj^i)N'N^ = jK, (25) 

which is used as the \ipdate equation for the boundary. 
The free surface moves if Eq. (25) does not hold for a 
given trial boundary. 

C. Second variation with respect to the surface 
parameterization 

In a first implementation we used a direct and explicit 
update algorithm moving the boundary into normal di- 
rection with a step-width that is determined by a param- 
eter T and the residual of Eq. (25). The discretization of 
this update can be found below in Eq. (55). Depending 
on the value of t this method exhibited strong instabili- 
ties as demonstrated below in Fig. 2. Although advanced 
techniques for determining an apt value for r seem to 
exist (cf. the program surface evolver by Brakke''^), we 
prefer a Newton-Raphson iterative method. This has the 
advantage of a faster convergence and a less strong de- 
pendence on T. A minor disadvantage is that it requires 
an additional variation of the surface free energy for the 
assembly of the geometry system. Using the same cal- 
culus as in Appendix B we find the second variation of 
the free energy contribution of a one-dimensional free 
surface, 

5'^F[St, St]=5(^J 7 gija°'Hli^6t]„ dA^ (26) 
= j^lSti^Qi^a'^^dti^dA (27) 

- / 7 {Kc 9^k a'^^ t%) {5t^^^ g^i a^'^ t\^) dA. 

J A 

For a two-dimensional surface the corresponding varia- 
tion contains two additional terms that are not given 

here for brevity. The last integral in Eq. (27) turns out 
to cause numeric instabilities in Newton's method. This 
is a rather surprising fact, because the calculation that 
led to Eq. (27) consists of two straightforward variations. 
If the last integral in Eq. (27) is omitted the algorithm 
becomes stable and accurate (cf. the tests in Sec. VI B). 

It is not only the free energy contribution F that de- 
pends on the shape of the surface. Also, the flow velocity, 
and by this, the viscous stress and the pressure depend 
on the shape. The formulation of the Newton method 
requires also the change of the fluidic stress integral due 
to changes of the free boundary, 

5(^1 aijN^fdA^ [6t]= J Seaij6N^[6t]dA 

+ [ SfaijN^ds/^[5t]di'+ [ 5f5aij[5i\N^ dA. (28) 

JE J A 
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The first two integrals on the right-hand side contain the 
changes of the normal vector (14) and the infinitesimal 
surface area dA = ^/adu due to changes of the bound- 
ary's shape. Both can be calculated along the lines of 
Appendix B. The third integral expresses the change of 
the fluidic stress (Ty at the boundary due to changes of 
its position. The shape changes are communicated to 
the flow and pressure fields via the boundary conditions 
(3) and (24) of the Stokes equations. Unfortunately, this 
very indirect response of the stress tensor on the changes 
of shape cannot be expressed exactly. We therefore have 
to assume that this term can expressed by derivatives of 
the stress tensor, i.e., 



5(Jij [6i\ » (Tij^kSt'^ 



(29) 



This means that the fluidic and the geometry system de- 
couple to the extent that the stress tensor in the vicin- 
ity of the boundary is not affected by small boundary 
changes. We note that this is not a consequence of split- 
ting the problem into two separate systems, but a general 
problem that equally applies to the combined approach. 
Altogether, the right-hand side of Eq. (28) becomes ap- 
proximately 



{SfaijN^){St';,gkia"%) 

- (^iV„tf„)a"''((5if^A^fe) + 6fa^,,kN'6t''} dA. (30) 

IV. DISCRETIZATION OF THE PROBLEM 

We implemented the above equations by means of a 
Galerkin approximation scheme which is known to work 
well for minimization problems. As variables we intro- 
duced the velocity components u and v in x- and y- 
direction, respectively, the pressure p, and additional 
variables r and s for the coordinates of the boundary 
parameterization vector t. The continuous fields are dis- 
cretized using ansatz functions, weighted with the corre- 
sponding degrees-of- freedom (DoF), 



d d 



(31) 



(32) 



where the sum runs over all DoFs. The fiuid velocity 
components u,v are approximated by the second-order 
finite elements (FEs) (/> and the pressure variable p by 
first-order FEs if). For the position variables r, s we have 
predominantly used second-order FEs, but for accuracy 
and other testing reasons we also tried first-order FEs. 
We denote the position FEs with x- AH FEs are of the 
Lagrange family,^ having ansatz functions that are 1 at 



exactly one node of the mesh and at all others. The 
DoFs are then equal to the function values at the nodes. 
This property is most convenient for the position vari- 
ables {rd, Sd) that coincide with the coordinates of the 
node d. 



A. The fluidic system 

The fluidic system is implemented in a standard way. 
Equation (2) is tested with the second-order FEs (p, while 
the continuity Eq. (1) is tested with the first-order FEs ip- 
We have implemented the following linear equation for 
the DoFs, which are collected to vectors u, v, p with com- 
ponents Ud,Vd,Pd respectively, 



K„ 



Kyy Kyp 





(34) 



with the entry matrices K and entry vectors L given by 

]de=vj '^(t>d-V(t)e dV - r]j> (t>d^-V(t)e dA, (35) 

V dV 

[Kvv]de =vj V.?id-V0e dV - ^d^-V^e dA, (36) 

V dV 

[Kup]de = - I {dxMi^edV + i MeN^dA, (37) 

JV JdV 



\i^vp\ de 



{dy(j)d)i^edV + f (l>dil'eNydA, (38) 

V JdV 



pu\de 



Ipddx^e dV, 





JV 


[^pv]de — 


- / -ipddyct^edV, 
Jv 


[Lu\d = 


/ (pdfxdV, 
Jv 


[Lv]d = 


Jv 



(39) 
(40) 
(41) 

(42) 



All integrals arc assembled in a loop over the elements 
and the sides of the mesh, using a fifth-order Gaussian 
quadrature rule. The fluidic system could likewise imple- 
ment the stationary Navier-Stokes equations with a small 
Reynolds number; we have chosen the Stokes equation for 
simplicity reasons here. 

The boundary conditions are imposed by a constraints 
technique for the matrix and for the right-hand side in 
Eq. (34). A constrained DoF Ud is expressed by an inho- 
mogeneity plus a weighted sum of other DoFs, 



Ud = Wd + '^ Wde 

e^d 



(43) 



which represent the boundary condition in question. The 
DoF Ud is then completely eliminated from the linear 



8 



system (34). By such constraint equations we imple- 
mented weak formulations of the kinematic boundary 
condition (3), i. e., 

= VKA^x+t^eiV,/) / MedA, (44) 
of the no-slip condition at the walls 

= y2{UeT^ + VeTy) [ (j>d(j)edA, (45) 

Jav 

and of the tangential projection of the free-surface bound- 
ary condition (24), 



in a discretized form as 



Tx N y -\-Ty N X 



2TyNy 



dV 



dA. (46) 



The constraint equations differ only in the values of Wde- 
The inhomogeneity Wd is zero in all three equations. 
Non-zero inhomogeneities would result, if also a surface- 
gradient term of the tension were taken into account in 
Eq. (24), or if the rigid walls performed a tangential 
movement. 

For the boundary condition (46), which is equivalent 
to an ideal slip condition, it is known that an improper 
choice of the normal direction can cause spurious contri- 
butions to the velocity field (see Behr,'^* Walkley et. al.,^^ 
and our remarks stated in the introduction). In the pres- 
ence of conservative forces only we did not find such spu- 
rious flows in our results. 

The fact that the formulation of the free boundary con- 
dition in terms of the DoF- constraints (46) cross-links 
all DoFs residing at boundary nodes, presents a serious 
problem. Each of the DoFs is in principle linked to all its 
neighbors on the boundary. This leads to a nearly filled 
system-matrix which is unfavorable regarding memory 
capacity and computing time. We found that an itera- 
tive method can overcome this problem. Instead of cross- 
linking a boundary DoF with all its neighbors, for some 
of them we take their old values, as is detailed in Ap- 
pendix D. After some iterations the full boundary con- 
dition (46) is established. The drawback of this scheme 
is that the constraint equations have to be re-assembled 
after every solution step of the fluidic system. 



B. The geometry system 

The geometry system employs a Newton method to 
perform the nonlinear search for the correct boundary 
position. This scheme corresponds to a minimization of 
the free energy F, while taking the fluidic stress into 
account. The boundary update equation can be written 



OF f 

= [Lr]dir, s,u,v,p) := — + / Xdcr'^' Nj dA, 
axd J A 

dF r 

= [Ls]d{f, s,u,v,p):= — + / XdCJ^'^Nj dA, 



(47) 



where d runs over the DoFs for each geometry variable, 
and L and F are understood as functions of the arrays 
r, s, etc. containing the DoFs. For the Newton-Raphson 
method the geometry system repeatedly has to solve the 
linear system of equations^" 



/ d[Lr\d/dre d[Lr\d/ds, ^ ^"'"^ 
V d[Ls]d/dre d[Ls]d/dse 



-*(ncw) 



.(old) 



-.(new) _ -.(old) 
•5e 



[Lr]d 

[Ls]d 



(old) 



(48) 



where r e [0, 1] is a step-size parameter. In all applica- 
tions we have used values of t between 0.1 and 1.0. 

The search for the correct boundary shape is strongly 
nonlinear in the position variables. In order to remove 
the main nonlinearities, which are caused by the surface 
metric expressions i/a and a"^ , the nodes of the elements 
are moved to their corresponding coordinates {rd, Sd) af- 
ter each step of the geometry system. Then, all integrals 
can be performed directly on the elements' edges. Also 
the normal vector can be taken from the elements' sides. 
In the previous section we used a convenient variational 
notation to express the change of the free energy contri- 
bution F by changes of the boundary parameterization. 
Essentially the same equations are obtained by differen- 
tiating the discrete version of F with respect to the DoFs 
which are the nodal degrees of freedom of the correspond- 
ing variables. The only difference is that the variation Sf^ 
in the continuous formulation must be replaced by the 
vectorial test-function Xd^ii and the variation 5t^^ by its 
tangential derivative Ta-Vx^e,. 



C. Controlling the tangential displacements of 
boundary nodes 

For a given discretization we must not only find the 
correct boundary shape, but its discretization should also 
remain well-proportionate. Very long and very short el- 
ement sides cause badly conditioned matrices and make 
the whole algorithm unstable. Several algorithms imple- 
menting the weak form of the free-surface boundary con- 
dition encounter these intrinsic instabilities of the bound- 
ary mesh. For the program surface evolver this mani- 
fests itself in shrinking and growing surface facets. It is 
therefore recommended to monitor the mesh quality and 
remove too small or split too large elements. Similar 
effects were reported by Brinkmann.^^ 

In Sec. IIIB the assignment of the boundary condi- 
tions to the fluidic and the geometry systems was de- 
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scribed. There, we found that the presence of incom- 
patible forces may easily destroy a free surface which es- 
sentially attempts to minimize the lengths j4(„) of the 
free-surface sides in each element m. Because all flu- 
idic stresses are constrained to have only normal com- 
ponents, we are free to use additional tangential force 
components for keeping the boundary mesh as regular 
as possible. This can be done during the assembly of 
the system matrices by weighting the surface tension by 
the element's side length divided by the average 

length of all element sides contributing to the free 

surface. Of course, this weighting factor becomes ineffec- 
tive if all sides have equal length. Any length difference 
of adjacent sides causes an additional force that tries 
to equalize them. The tension forces for each element 
side are then equivalent to a first variation of the func- 
tional 7A^^^/(2(A(m))), which describes a rubber band 
with Hookean forces. Instead of (5F[(5t] from Eq. (17) we 
thus assemble on each element 



2(^(„,)) 



M(„)[(5t]. 



(49) 



The second variations of and are not pro- 

portional to each other, 



2(^(m)) 



A 



1' 



(m) 



{5A^^)[5t]f (50) 



In the implementation we therefore took only the first 
term on the right-hand side of Eq. (50). In this sense we 
did not strictly implement the behavior of a rubber band, 
but yet a stabilized version of the free-surface tension 



terms. After convergence all boundary sides of the mesh 
representing the free surface have equal lengths and the 
extra terms Ai^^n) / {^{m)) do not change the behavior of 
the free surface. 



V. SUMMARY OF THE ALGORITHM 

Here, wc provide a short overview of the complete al- 
gorithm. The required steps are as follows: 

1. Choose an initial mesh and initial ambient 
pressure po. 

2. Until convergence repeat the following steps: 

(a) Smooth the inner mesh if it is too distorted. 

(6) Repeatedly solve the fluidic system for p, 
u and V, until the slip boundary condition is 
established. 

(c) Subtract the average from p. 

(d) Solve the geometry system for the new 
boundary. At the same time search for the 
value of pq that keeps the volume unchanged. 

(e) Set the mesh boundary nodes to the 
parameterization values of the geometry 
system. 



The fluidic system is assembled according to Eqs. (34)- 
(42) with constraints that account for the proper bound- 
ary conditions. To give the full algorithm at this point, 
we summarize also the terms of the geometry system. 



The update Eq. (48) is written as 



Ksr K,, I \ s (-W) ) - ^ [ lJ ^ [ [ s (old) 



with entries that are assembled per element m, 



(Am)} 



(51) 



(52) 



[A-(;")]de = - / XdXe(e.-Vp)(e,.N) dA+ J Xd (Vxe-T){(e,.aN)(e,.T) - (e,.aT)(e,.N)} dA (53) 

^("») ^(m) 'v4, s r 

+ 7a^ i^Xd-T){Vxe-T)dA 

[K(^^^^ = - J XdXe{e.-Vp){ey-N) dA+ J Xd (Vxe-T){(e,.aN)(e,.T) - (e,.<7T)(e^.N)} dA. (54) 



The remaining entries can be obtained by permutations of x and y together with r and s. Again, constraints have 
been used to keep the contact-lines pinned. 
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VI. NUMERICAL EXPERIMENTS 

We performed all our test cases for a two-dimensional 
fluid. The programs were written using the open-source 
C++ library libmesh'^^ which allows to change the ele- 
ments' geometry in a user's routine and provides a pow- 
erful constraint method. 



A. The instability of a "direct explicit update" 
algorithm 

In our first numerical example wc do not use Newton's 
method with the update-rule (51), but instead with the 
direct and explicit update 




The stability of this update rule delicately depends on the 

step-size parameter r. The allowed range of t strongly 
depends on the size of the elements, the curvature, etc. 
Figure 2 depicts the most simple situation where a ho- 
mogeneous prcssiire field deforms the boundary into a 
circular arc with radius R = —1/k — Poll- The update- 
rule (55) is stable for 14 FEs and a given step-size while 
it is unstable for the same step-size with 52 FEs. 



B. Testing the accuracy of the Newton algorithm 

In order to confirm the accuracy of the curvature ap- 
proximation we have tested two cases that can be solved 
analytically. Similar to the calculation in Appendix A, a 
prescribed pressure determines the free surface's shape. 
Then, the approximation in Eq. (29) becomes exact and 
simplifies to 

5cFij = -gijP,kSt''. (56) 

Thus, all possible approximation errors must be due to 
the discretization of the curvature. 

Figure 3a depicts the most simple situation where a 
homogeneous pressure field deforms the boundary into a 
circular arc, as in the previous example. The free surface 
shape is approximated by the sides of 5 second-order FEs. 
In dimensionless units the surface tension is 7 = 1, and 
the prescribed pressure po = 2 produces as the exact so- 
lution a circle with radius R = 1/2. The relative error of 
the numerically resulting radius, and thus also of the cur- 
vature, is about 8.6 x 10~^. This value has been obtained 
from the position of the topmost node. An alternative 
approach for calculating the approximation error is vi- 
sualized in Fig. 3b. We calculated the normal vectors at 
each node from the resulting finite-element's side. Due to 
the elements being second-order we got a single normal 
vector for second-order nodes. At vertices, where two ele- 
ments meet and where the surface parameterization is not 




-0.4 -0.2 0.2 0.4 

X 




-0.4 -0.2 0.2 0.4 



X 

FIG. 2: The stability of a direct explicit update algoritlim 
strongly depends on the ratio of the step-size r and the cl- 
ement size. We have used the same t = 0.05 for two dif- 
ferent numbers of approximating FEs (first-order). A pre- 
scribed homogeneous pressure po = 2 is applied which bends 
the free surface into a half-circle with radius 1/2 (with 7=1). 
Panel (a) depicts the converged result for 14 FEs after more 
than 500 steps. Panel (b) siiows a mesh with 52 FEs after 
only 12 steps. For this combination of step-size and element 
size the direct update algorithm is unstable, and the mesh was 
completely destroyed after a few more steps. Using second- 
order FEs the instability was similar. In both panels the 
dashed half-circles indicate both, the exact solution and the 
initial geometry. 



smooth, we averaged the two normal vectors. The cur- 
vature estimate at a node in Fig. 3b is then given by the 
curvature radius of a circle that connects the two neigh- 
bors of the specific node, given their appropriate normal 
vector. Thus, we explicitly reconstructed the curvature 
from the change of the normal vector along the surface. 
It is clear by construction that the normal vector of the 
contact-nodes cannot be correctly estimated. This causes 
the four outliers in Fig. 3b. All other nodes fit well. 

A comparison with Fig. 2a, where the result of a first- 
order approximation can be seen, makes clear that it is 
crucial to use a second-order parameterization. The rel- 
ative error of the curvature in Fig. 2a is 5.0 X 10~^, three 
magnitudes larger than in Fig. 3a. 

In the next accuracy test, depicted in Fig. 4, the pres- 
sure is still prescribed, but it varies in space. As above, 
we apply a pressure for which the resulting boundary 
shape is known. Figure 4a illustrates the approximation 

of a sinusoidal boundary height function y = h{x) = 
asm{Px) that is caused by the corresponding pressure 
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FIG. 3: A prescribed homogeneous pressure po = 2 bends the 
free surface into a half-circle with radius 1/2 (with 7 = 1). 
Panel (a) presents the approximation with only five second- 
order FEs. The boundary nodes are indicated by circles (ev- 
ery second is a second-order node). The exact solution is 
indicated by the dashed half-circle, while the starting geome- 
try was the straight connection between the fixed endpoints. 
Good convergence was reached after 100 iterations with a 
step-size parameter t — 1. The topmost node misses its exact 
position only by a relative error of only 8.6 x 10~^ . This is also 
the error of the overall curvature approximation. In panel (b) 
an estimate of the curvature k was obtained by reconstructing 
the normal vectors and their change along the surface from 
the boundary shape of the FEs. This estimate compares very 
well with its expected value of —2.0. The outliers near the 
contact points are artefacts due to the reconstruction of the 
normal vectors. 



field 



a0^ sin{(3x) 



7[l + a2/32 cos2(/3a;)]3/2' 



(57) 



Again, the approximation in Eq. (29) becomes exact, and 
we expect the same discretization errors as in the previ- 
ous example. The curvature's relative error is larger than 
in the previous example because the curvature is bigger 
compared to the number of nodes. Nevertheless, the er- 
ror is still small enough to return the expected boundary 
shape within reasonable accuracy. It decreases with the 
number of approximating elements. If a first-order ap- 
proximation is used it is much larger, maybe intolerably 
large. 

Concerning the discretization errors of the curvature 
the accuracy test in Fig. 4 covers already the general case. 
According to the construction of the algorithm the flow 




0.4 0.6 



FIG. 4: In panel (a) the expected sinusoidal boundary shape 
y = h(x) = 0.25 sin(47ra::) with surface tension 7 = 1 is well 
recovered by 40 second-order FEs. The shape is generated by 
the prescribed pressure of Eq. (57) . As in Fig. 3 the nodes are 
indicated by circles, the exact solution by the dashed curve 
and the initial geometry was the straight connection between 
the fixed endpoints of the surface. Good convergence was 
reached after 60 iterations with a step-size parameter t = 1. 
The nodes' position at the maxima is off by a relative er- 
ror of 6.7 X 10"^ (2.3 X 10"^ for 80 FEs and 2.5 x 10"* for 
120 FEs, not shown). The elements' side-lengths vary only by 
±0.007%. This small deviation demonstrates that the mesh 
regularization method does not influence the flnal behavior 
of the free boundary. The approximation quality is dramati- 
cally worse for 40 first-order FEs, yielding an estimated error 
of 2.0 X 10"^ (not shown). In Panel (b) the applied pressure 
is compared with curvature estimated by a reconstruction of 
the normal vectors. The large deviations at the extrema do 
not affect the overall approximation of the sinusoidal shape. 



exerts stress on the boundary only in normal direction. 
It makes no difference whether this stress is of viscous 
nature or due to a pressure difference. 



C. A deformed micro-droplet 

In order to explicitly show that the quality of the cur- 
vature discretization does not depend on the origin of 
the applied normal stress we like to return to the intro- 
ductory motivation for the present work. The previous 
examples were analytically solvable. The form and in- 
ternal streaming of micro-droplets, however, cannot be 
determined analytically. 
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Non-conservative force density: field lines 
and magnitude (in N/m^) 




Conservative force potential (in N/m^) 




FIG. 5: The force density that models the effect of the SAW 
in the droplet given in Fig. 6. Panel (a) depicts the non- 
conservative part that causes the flow; (b) shows the potential 
of the conservative part that contributes only to the pressure. 
The same non-conservative force density has also been used 
in Fig. 7. 



In the experiment the internal flow is agitated by a 
surface-acoustic wave (SAW) due to the acoustic stream- 
ing effect.^"^ Because the very details of the SAW's impact 
are not known, we here model it by a body force that is 
active in the fluid only, as depicted in Fig. 5. The force 
is concentrated in a narrow channel starting at the left 
contact point where the SAW hits the fluid and contin- 
uing into the fluid. It essentially carries the fluid along 
this channel, from the entry point of the SAW into the 
droplet, giving rise also to a back-flow. ^■^ Additionally, 
the force has a strong conservative portion that is bal- 
anced by the pressure in the fluid. 

The resulting stationary droplet shape and the inter- 
nal velocity and pressure fields are presented in Fig. 6. 
The initial shape was a half-circle with the same two- 
dimensional volume. The deformed boundary consists 
of two regions, one with negative curvature (as the ini- 
tial half circle) and another one at the right flank of the 
droplet with positive curvature. The material proper- 
ties are those of water and air at room temperature, i. e. 
7j = IQ-^kg/ms and 7 = 72.8 x 10"^N/m. The - admit- 
tedly strange - deformation of the droplet qualitatively 
agrees with the experimentally observed jumping droplet 
in Fig. 4 of the publication by Wixforth et. al.}^ The de- 
formation is due to the large conservative contribution of 
the driving force and the resulting pressure. The viscous 
forces for the given velocities are far too weak to lead to a 
substantial deformation of the free surface. The capillary 
number for the illustrated flow is Ca k, 10~^, the Bond 




Velocity streamlines 
and magnitude (in 10~*m/s) 
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FIG. 6: A deformed micro-droplet, sitting on a flat substrate 
with pinned contact points. The deformation is due to an 
internal pressure and viscous flow, both caused by the body 
force density illustrated in Fig. 5. The material properties 
are those of water surrounded by air at room temperature. 
Its two-dimensional "volume" is that of the initial half-circle 
with radius 0.28 mm. Panel (a) illustrates the computational 
grid, consisting of second-order elements. The side-lengths 
of the free-surface facets differ only by 4.5 x 10~^%. Panels 
(b) and (c) depict the flow and the pressure, respectively. 
Note that the deformation is predominantly caused by the 
pressure which corresponds to the case that Ca <C Bo. Good 
convergence was reached after 7 iterations with a step-size 
parameter r = 0.5. 
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Velocity streamlines 
and magnitude (in 10~''m/s) 




1 1 1 1 1 1 1 1 1 1 1 1 

^ <f ' ■ 


■ + T-o-N 
- X N-o-N 

1 


■ o 

V * ^>^^ 

** - 

1,1, ^ ,1,1 



0.004 
0.002 


-0.002 
-0.004 
-0.006 

10 20 30 40 50 60 

Free-surface node number 

FIG. 7: A similar micro-droplet as in Fig. 6, deformed only 
by the viscous stress at the boundary. The flow is driven by 
the non-conservative force density depicted in Fig. 5a. The 
conservative part of the force vanishes such that the pressure 
is constant inside the droplet. In order to obtain a compa- 
rable deformation we have used a 10^ times smaller surface 
tension than that of a water-air interface. This corresponds 
to the case Bo <C Ca. Good convergence was achieved after 
30 iterations with r = 0.1. 



number is around one. Although the free surface is sig- 
nificantly deformed, its discretization by finite-element 
sides is as regular as possible. Their lengths vary only 
by 4.5 X 10~^%. This guarantees that the behavior of 
the boundary is indeed that of a free surface and is not 
disturbed by the automatic regularization technique de- 
scribed in Sec. IV C. Figure 6d quantifies the normal 
stress condition. For each node we integrated the nor- 
mal and the tangential component of the normal stress, 
weighted with the corresponding ansatz function of the 
node. The tangential component vanishes perfectly. The 
normal component coincides well with the reconstruction 
estimate of the curvature as in the previous examples. 
Thus, the free-surface boundary condition is indeed sat- 
isfied. 



In order to prove that our algorithm can likewise pro- 
duce stable results in the parameter regime Bo <^ Ca 
we consider a droplet that is deformed only by viscous 
stress at the boundary. In Fig. 7 we have used the same 
non-conservative force that is visualized in Fig. 5, but 
we omitted the conservative part. Thus, the pressure 
was constant and Bo = 0. With the same water-air in- 
terface tension as in the previous example the droplet 
would hardly be deformed. To obtain a comparable de- 
formation as in the previous case together with Ca « 1 
we took an artificial 10^ times smaller surface tension. 
In this example the stress that deforms the free surface 
depends much stronger on the shape of the surface itself. 
Thus, the approximation in Eq. (29) becomes question- 
able. It was necessary to reduce the step-size parameter 
T to a smaller value than in the previous examples. 



VII. SUMMARY AND OUTLOOK 

Within this work we presented a weak formulation of 
free-surface boundary problems in arbitrary coordinate 
systems. The steps of the derivation are physically and 
mathematically founded using variational techniques for 
the Stokes equations and the differential geometry of 
the surface. We found that the applicability of differ- 
ent numerical treatments for the curvature terms depend 
strongly on the scales of the system. Our method is de- 
signed for Bond and capillary numbers assuming values 
from zero up to unity. Which one is larger plays no role. 

A decisive benefit of our method is the automatic con- 
trol of mesh regularity at a free surface. Many algorithms 
implementing the weak form of the free-surface boundary 
condition encounter intrinsic instabilities of the bound- 
ary mesh. Often, it is therefore necessary to create a 
completely new mesh after several iteration steps. Our 
formulation includes a smooth transition to the behav- 
ior of a rubber band when the boundary mesh becomes 
distorted. This leads to an inherent regularization of the 
mesh without affecting the behavior of the free surface. 

As another important result we find that for physical 
reasons the geometry variables for the parameterization 
of the free surface should be approximated on the same 
level of accuracy as the velocity variables. This substan- 
tiates numerical observations reported by Bansch.^^ 

The quality of our numerical approach is tested by two 
analytically solvable examples. We explicitly plot the 
curvature of the free surface and the stress that causes the 
deformation. This confirms that the free surface bound- 
ary condition is indeed satisfied, not only in the weak 
sense which is implemented but even leads to a reliable 
reconstruction of the curvature by the normal vectors of 
the finite element's sides. Two further examples illustrate 
that the ratio of capillary number and Bond number has 
only a weak influence on the stability of the algorithm. 

The presented covariant formulation opens the possi- 
bility to utilize the powerful differential geometric de- 
scription of free surfaces in finite-element implementa- 
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tions of the Stokes equations. It thus provides a natural 
approach to treat surfaces and interfaces with a richer 
behavior such as lipid vesicles containing bending stiff- 
ness, area constraints, and much more. Many potential 
applications can be found in the hterature on lipid vesi- 
cle geometry, where more complicated expressions for the 
surface's free energy contribution are in use."*^"^^ 

Extensions of the presented approach towards moving 
contact-lines, towards time-dependent flows and a three- 
dimensional implementation are possible. There are still 
some hurdles to be overcome that can be clearly seen in 
our derivation. One of them is the principally unknown 
mutual dependence of the stress tensor and the surface 
parameterization where we had to introduce the approx- 
imation (29). Another one is the imderstanding of the 
numeric instabilities caused by the last integral in sec- 
ond variation of the surface's free energy (see Eq. (27)). 

These extensions would also provide a solid basis for 
the theoretical understanding of particle transport in 
surface-acoustic- wave-driven flows. ^^'^^"^^ 
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The Euler-Lagrange equation for finding the extremal F 
by varying h is then 

dJ" d dJ^ d dT 



° dh dxd{d^h) dyd{dyh) ^^^^ 
= <^>{x,y,h{x,y))--fK{x,y) (A6) 

where k is the curvature of A, given by 
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Because the pressure is given by the potential, the 
Euler-Lagrange equation is equivalent to the free-surface 
boundary condition for a static fluid. 



-p{x, y, h{x, y))+po= jk{x, y). 



(A8) 



APPENDIX A: A STATIC DROPLET IN 
CARTESIAN COORDINATES 

The aim of this appendix is to recall the variational 

techniques in a simple three-dimensional Cartesian setup 
before going to arbitrary coordinates in appendix B. The 
argumentation is similar to the one given by Cuvelier.^^ 
We describe the special case where the fluid's two- 
dimensional free surface A can be described by a height 
function 



At this point it is easy to see that pq plays the role of a 
Lagrange multiplier for a volume constraint. Adding the 
term 



W = A / dxdy h{x, y) 
Je 



(A9) 



to F gives an additional constant A in the Euler-Lagrange 
equation, just as the pressure offset po- Because po is yet 
undetermined we may identiiy it with A. 



A: z = h{x, y) 



(Al) 



which is non-zero over a certain region {x,y) G E. The 
Stokes equations for the static situation with a conserva- 
tive force fi = — reduce to 







-P,i - ^, 



(A2) 



with the solution p(x) = po — ^{x). The undeter- 
mined homogeneous term po will be identified as the La- 
grange multiplier for the constraint of constant volume 
V = fj^dxdyh{x,y). 

The free energy of the system consists of the surface- 
integral of the constant surface tension and the volume- 



APPENDIX B: VARIATIONAL CALCULUS FOR 
THE SURFACE'S PARAMETERIZATION 

In order to prove equality (17) we express the change 
of the surface's free energy functional (9) by the change 
of the Jacobi determinant a/o of the surface parameteri- 
zation. With the infinitesimal surface area dA = ^/adu 
the variation of the surface free energy becomes 



5F[5t] = s(^J jdA^ ^ J 7 S^dv 



Ie^ dt\, 



5t\ du. (Bl) 
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The dependence of ^Ja on the tangent vectors follows 
from its definition as the determinant of the surface met- 
ric. For a two-dimensional surface it reads 



an oi2 

0-1\ 0,22 



1 



= :^e"'^^'''gijgut:^t>^0t%tls (B2) 



where e"^ is the permutation symbol in two dimensions, 
= 4-1 = 1 R = 2 (B3) 



a = (3 
+1 a = 1, 13 
-1 a = 2, (3=1 



which is a relative surface tensor with weight +1. The 
absolute tensor results as 



-a/3 



(B4) 



This is analogous to the completely antisymmetric tensor 

in three dimensions, described in detail by Aris.^^ With 
the antisymmetric tensor we obtain the inverse surface 
metric as 



A formal derivative of (B2) yields 
da 

1 da 



dy/a 

dt\ ~ 2y/adt\ 



,/3 



(B5) 

and (B6) 
(B7) 



which can be inserted into (Bl) to give the desired re- 
sult (17). 

For a one-dimensional curve in two-dimensional space 
the same formula can be derived, but the notation may be 
somewhat confusing. Summation over the single surface 
index makes no sense, nevertheless, wc still have to dis- 
tinguish between co- and contravariant relative tensors, 
i. e., 

a = On = gijt\^t\i (B8) 
a^^ = 1/aii because a^^an = a°'^actp = 1. (B9) 
The formal derivative then becomes 



^^/a 
'di' 



1 



1 



- = 7r7=^9ijt\ = \fagijt\ — = y/agijt\a^^ 
1 ^v'* Oil 

(BIO) 

which completes the result for the one- dimensional sur- 
face. 



APPENDIX C: INTEGRATION BY PARTS OF 
THE TENSION FORCES 

In order to see that Eq. (18) follows from Eq. (17) we 
remove the surface covariant derivative from Sf^ by an 



integration by parts and obtain 



6F[6t] = - / ja''H'„^gij6t^ - / j^f.a^f'f^g.jSt^ 

J A J A 

+ i jiy^a^^f^gijSt^ (CI) 

JdA 

where the covariant surface vector i^p is tangential 
to A and normal to dA. We can express the surface- 
derivatives t^^^p by the tensor b^p of the second funda- 
mental form of the surface from Eq. (16) (cf. to Aris,^^ 
p. 216), 



arriving at 



^boc^N' = kN\ 



(C2) 



(C3) 



Wc have used the definition of the curvature as the 
trace of the tensor of the second fundamental form as in 
Eq. (15). For a two-dimensional surface this is twice the 
mean curvature n = 2H = a^^ba/s, for a one-dimensional 
surface we have only one entry k = a^^bn. 

As consistent with the standard literature, ^^'^^ the 
term 5F from Eq. (17) comprises a curvature term in 
normal direction 



(C4) 



and a term accounting for the surface-gradient of 7. The 
space vector 



7,/3 



(C5) 



is tangential to the surface. The third term on the right- 
hand side of Eq. (CI), which is an integral over the 
contact-line dA, vanishes because for a pinned droplet 
6t^ = vanishes on the contact-line. 



APPENDIX D: INVOKING CONSTRAINTS FOR 
THE SLIP BOUNDARY CONDITION 

The tangential components of the free-surface bound- 
ary condition correspond to a slip boundary condition. 
When this condition is expressed as a set of constraints 
for the DoFs, we obtain one equation like (46) per each 
DoF at the free surface. Because the derivatives of 
the ansatz functions (j) from (31) generally do not van- 
ish at proximate nodes, the constraint equations contain 
non-vanishing weights for all DoFs that arc located on 
the same element. Therefore, the constraints for DoFs 
that are connected to two adjacent elements create inter- 
dependencies of DoFs also on other elements. This is 
illustrated in Fig. 8a. As a result, all DoFs on the free 
surface implicitly depend on each other. After the con- 
straints are re-sorted such that DoFs are constrained only 
in terms of non-constrained ones, it turns out that the 
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FIG. 8: A sketch of the cross-dependencies among the DoFs 
located on three elements. The free surface is indicated by 
the thick curve. Constrained DoFs are surrounded by small 
circles. The nodes carrying the corresponding constraining 
DoFs are surrounded by curves drawn in the same style (solid 
and dashed for vertices; dotted and dash-dotted for second- 
order nodes). Panel (a) depicts the full inter-dependencies 
while in (b) the DoFs located at vertices depend only on DoFs 
located at inner nodes. By taking the values of the missing 
adjacent DoFs located on the free surface as inhomogeneities 
instead of constraints in (b) constraints are decoupled. The 
correct constraint equations are then established after some 
iteration steps. 

free-surface DoFs depend on all DoFs in the element layer 
near the surface. 



As a strategy to avoid this full dependency we replace 
the constraint equation of type (43) by 

Ud^Wd+ ^ WdeUe + ^ Wdeu'f^'^\ (Dl) 

where the sums run over two complementary sets A^; and 
Ad- The DoFs in A^; contribute to the constraint for Ud 
in the usual way, while those in A^ have been substituted 
by their old values ui"'"^"* and thus contribute to the in- 
homogeneity. There is some freedom in the choice, which 
of the participating DoFs in one element are in A^ and 
which are taken into A^. We found that the combination 
illustrated in Fig. 8b works well: For the DoFs located at 
element vertices we take the DoFs that belong to adjacent 
nodes on the free surface as inhomogeneities; all other 
constraining DoFs are located at inner nodes and are not 
constrained. The DoFs located at the second-order nodes 
on the free surface acquire their full constraints. When 
all constrained DoFs are expressed by non-constrained 
DoFs, then the resulting constraint equations will only 
contain DoFs that are located at inner nodes of three 
adjacent elements. This presents a sufficient decoupling 
of the constraint equations to yield an efhcient algorithm. 

Although the boundary condition given by Eq. (Dl) 
is not the correct one when the true velocity field has 
not yet been determined, it still improves as the velocity 
field tends to the proper solution. Thus, there is hope 
that the correct boundary condition is established by the 
successive use of Eq. (Dl) using increasingly good val- 
ues for the values In numerical experiments the 
scheme for splitting the cross-dependencies as illustrated 
in Fig. 8b turned out to be the only one that works. 
In the examples of Figs. 6 and 7, it took about 20 iter- 
ation steps to establish the correct boundary condition 
from scratch, and 5 iteration steps to re-establish it after 
a change of the mesh. This could be readily observed 
because after the first iteration step the velocity field ex- 
hibited oscillations at the boundary nodes that ceased 
during iteration. 
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